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Magnetic turbulence in compressible fluids 

By Jungyeon Cho f and A. Lazarian $ 



We present high resolution numerical simulations of compressible magnetohydrody- 
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namic (MHD) turbulence. We concentrate on studies of spectra and anisotropy of veloc- 
ity and density. We describe a technique of separating different magnetohydrodynamic 
(MHD) modes (slow, fast and Alfven) and apply it to our simulations. For the mildly su- 
personic case, velocity and density show anisotropy. However, for the highly supersonic 
case, we observe steepening of the slow mode velocity spectrum and isotropization of 
density Our present studies show that the spectrum of density gets substantially flatter 
than the Kolmogorov one. 



1. Introduction 

Most astrophysical systems, e.g. accretion disks, stellar winds, the interstellar medium 
(ISM) and intercluster medium are turbulent with an embedded magnetic field that in- 
fluences almost all of their properties. This turbulence which spans from km to many kpc 
(see discussion in Armstrong, Rickett, & Spangler 1995) holds the key to many astro- 
physical processes (e.g., transport of mass and angular momentum, star formation, frag- 
mentation of molecular clouds, heat and cosmic ray transport, magnetic reconnection). 
Q-f Statistics of turbulence is also essential for the cosmic microwave background (CMB) 

radiation foreground studies (Lazarian & Prunet 2002). In this brief, using high resolu- 
tion simulations, we discuss statistics of 3D MHD turbulence and present new results on 
density fluctuations. 

Why do we expect astrophysical fluids to be turbulent? A fluid of viscosity v becomes 
turbulent when the rate of viscous dissipation, which is ~ v/L 2 at the energy injection 
scale L, is much smaller than the energy transfer rate ~ Vl/L, where Vl is the velocity 
dispersion at the scale L. The ratio of the two rates is the Reynolds number Re = VlL/v. 
In general, when Re is larger than 10 — 100 the system becomes turbulent. Chaotic 
structures develop gradually as Re increases, and those with Re ~ 10 3 are appreciably 
less chaotic than those with Re ~ 10 8 . Observed features such as star forming clouds are 
very chaotic for Re > 10 s . This makes it difficult to simulate realistic turbulence. The 
currently available 3D simulations containing 512 grid cells along each side can have Re 
up to ~ O(10 3 ) and are limited by their grid sizes. Therefore, it is essential to find "scaling 
laws" in order to extrapolate numerical calculations (Re ~ O(10 3 )) to real astrophysical 
fluids (Re > 10 8 ). We show below that even with its limited resolution, numerics is a 
great tool for testing scaling laws. 

Kolmogorov theory provides a scaling law for incompressible non-magnetized hydro- 
dynamic turbulence. This law provides a statistical relation between the relative velocity 
vi of fluid elements and their separation I, namely, vi ~ Z 1 / 3 . An equivalent description is 
to express spectrum E(k) as a function of wave number k (~ 1/0- The two descriptions 
are related by kE(k) ~ vf. The famous Kolmogorov spectrum is E(k) ~ fc~ 5 / 3 . The 
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applications of Kolmogorov theory range from engineering research to meteorology (see 
Monin & Yaglom 1975) but its astrophysical applications are poorly justified and the 
application of the Kolmogorov theory can lead to erroneous conclusions (see reviews by 
Lazarian et al. 2003; Lazarian & Yan 2003) 

Let us consider incompressible MHD turbulence first. There have long been an under- 
standing that the MHD turbulence is anisotropic (e.g. Shebalin et al. 1983). Substantial 
progress has been achieved recently by Goldreich & Sridhar (1995; hereafter GS95), who 
made an ingenious prediction regarding relative motions parallel and perpendicular to 
magnetic field B for incompressible MHD turbulence. An important observation that 
leads to understanding of the GS95 scaling is that magnetic field cannot prevent mixing 
motions of magnetic field lines if the motions are perpendicular to the magnetic field. 
Those motions will cause, however, waves that will propagate along magnetic field lines. 
If that is the case, the time scale of the wave-like motions along the field, i.e. ~ 1\\/Va, 
is the characteristic size of the perturbation along the magnetic field and Va = B/^JAirp 
is the local Alfven speed) will be equal to the hydrodynamic time-scale, l±/vi, where 
Zj_ is the characteristic size of the perturbation perpendicular to the magnetic field. The 

1/3 

mixing motions are hydrodynamic- likef. They obey Kolmogorov scaling, vi oc ij_ , be- 
cause incompressible turbulence is assumed. Combining the two relations above we can 
get the GS95 anisotropy, In oc (or fcii oc in terms of wave- numbers) . If we inter- 
pret l\\ as the eddy size in the direction of the local magnetic field and l± as that in the 
perpendicular directions, the relation implies that smaller eddies are more elongated. 

GS95 predictions have been confirmed numerically (Cho & Vishniac 2000; Maron & 
Goldreich 2001; Cho, Lazarian & Vishniac 2002a, hereafter CLV02a; see also CLV03); 
they are in good agreement with observed and inferred astrophysical spectra (see CLV03). 
However, the GS95 model considered incompressible MHD, but the media in molecular 
clouds is highly compressible. Does any part of GS95 model survives? Literature on the 
properties of compressible MHD is very rich (see reviews by Pouquet 1999; Cho & Lazar- 
ian 2004 and references therein). Higdon (1984) theoretically studied density fluctuations 
in the interstellar MHD turbulence. Matthaeus & Brown (1988) studied nearly incom- 
pressible MHD at low Mach number and Zank & Matthaeus (1993) extended it. In an 
important paper Matthaeus et al. (1996) numerically explored anisotropy of compressible 
MHD turbulence. However, those papers do not provide universal scalings of the GS95 
type. After the GS95 model, Lithwick & Goldreich studied scaling relations for high-&e£a 
plasmas and Cho & Lazarian (2002; hereafter CL02) for \ow-beta plasmas. 

The complexity of the compressible magnetized turbulence with magnetic field made 
some researchers believe that the phenomenon is too complex to expect any universal 
scalings for molecular cloud research. High coupling of compressible and incompressible 
motions is often quoted to justify this point of view. 

Below we shall provide arguments that are suggestive that the fundamentals of com- 
pressible MHD can be understood and successfully applied to astrophysical fluids. In 
many astrophysical fluids the regular (or uniform) magnetic field is comparable with 
the fluctuating one. Therefore for most part of our discussion, we shall discuss results 
obtained for SV <~ dB/y/Airp ~ Bo/y^Awp, where 5B is the r.m.s. strength of the random 
magnetic field. 

Our work done during the summer program at CTR that we report here was an 
attempt to get a better insight into the physics of MHD turbulence using high resolution 

f Recent simulations (Cho et al. 2003) suggest that perpendicular mixing is indeed efficient 
for mean magnetic fields of up to the equipartition value. 
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simulations (i.e. simulations within a 512 3 box). In what follows we discuss our numerical 
approach (sect. 2), our results obtained for the velocity and magnetic field statistics (sect. 
3) and the density statistics (sect. 4). We concentrate on studies of spectra and turbulence 
anisotropy. 



2. Numerical approach 

We use a third-order accurate hybrid essentially non-oscillatory (ENO) scheme (see 
CL02) to solve the ideal isothermal MHD equations in a periodic box: 

dp/dt + V • (pv) = 0, (2.1) 
dv/dt + v Vv + p- 1 V(a 2 p) - (V x B) x B/Anp = f, (2.2) 

dB/dt - Vx(vxB)=0, (2.3) 

with V • B = and an isothermal equation of state. Here f is a random large-scale 
driving force, p is density, v is the velocity, and B is magnetic field. The rms velocity 
5V is maintained to be approximately unity, so that v can be viewed as the velocity 
measured in units of the r.m.s. velocity of the system and B / y/4irp as the Alfven velocity 
in the same units. The time t is in units of the large eddy turnover time (<~ L/SV) and 
the length in units of L, the scale of the energy injection. The magnetic field consists of 
the uniform background field and a fluctuating field: B = Bo + b. 

We drive turbulence solenoidally in Fourier space and use 512 3 points, Va = Bo/y/^np = 
1, and po = 1. The average rms velocity in statistically stationary state, SV, is around 1. 

For our calculations we assume that B /\/4jrp ~ 5B/\JAirp <~ SV. In this case, the 
sound speed is the controlling parameter and basically two regimes can exist: supersonic 
and subsonic. Note that supersonic means low-beta and subsonic means high-beta. When 
supersonic, we consider mildly supersonic (or, mildly low-/3) and highly supersonic (or, 
very low-/3)f . 

2.1. Separation of MHD modes 

Three types of waves exist (Alfven, slow and fast) in compressible magnetized plasma. 
The slow, fast, and Alfven bases that denote the direction of displacement vectors for 
each mode are given by 

is oc (-1 + a - VD)k\\£\\ + (1 + a - VD)k±k±, (2.4) 

if oc (-1 + a + \/^)fc||k|| + (1 + a + VD)k ± k ± , (2.5) 

i A = -<P = 'k± x k,|, (2.6) 

where D = (1 + a) 2 — 4acos#, a = a 2 jV\ = (3(j/2), 9 is the angle between k and B , and 
(p is the azimuthal basis in the spherical polar coordinate system. These are equivalent 
to the expression in CL02: 



£ s oc fc||ki| + 



l-VD-P/2 [fc|| 
1 + VD + (3/2 [kj_ 



2 



fciki, (2.7) 



f The terms "mildly" and "highly" are somewhat arbitrary terms. We consider these two 
supersonic cases to cover a broad range of parameter space. Note that Boldyrev, Nordlund, & 
Padoan (2002) recently provided a Mach number dependence study of the compressible MHD 
turbulence statistics where only two regimes are manifest: essentially incompressible and essen- 
tially compressible shock-dominated (with smooth transition at some M 3 of order unity). 
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Figure 1. Separation method. We separate Alfven, slow, and fast modes in Fourier space by 
projecting the velocity Fourier component Vk onto bases £a, £ s , and £/, respectively. Note that 
£a = —0. Slow basis £ s and fast basis lie in the plane defined by B and k. Slow basis £ s lies 
between —6 and ky. Fast basis £/ lies between k and k±. From Cho & Lazarian (2003). 



l-y/D + 0/2 
SJ 1 + V^D - (3/2 

(Note that 7 = 1 for isothermal case.) 

Slow and fast velocity components can be obtained by projecting velocity Fourier 
component Vk onto £ s and £/, respectively. See Cho & Lazarian (2003; hereinafter CL03) 
for discussion regarding how to separate slow and fast magnetic modes. We obtain energy 
spectra using this projection method. When we calculate the structure functions (e.g. for 
Alfven modes) in real space, we first obtain the Fourier components using the projection 
and, then, we obtain the real space values by performing Fourier transform. 

We tested our technique of separation in CL03 for a case when the separation is possible 
in real space and got essentially identical results with our statistical technique. Therefore 
we believe that our separation procedure works reliably. 



- z 



fciikii + k±k A 



(2.8) 



3. Velocity and magnetic field spectra 

We show in Fig. 2 new results from high resolution simulation of highly supersonic 
magnetically dominated media. The sonic Mach number, M s , is ~ 10. Fig. 2 shows that 
most of the scaling relations we previously found in the low resolution simulations are 
still valid in the high resolution simulation. Especially anisotropy of Alfven, slow, and 
fast modes is almost identical to the one in the previous studies. However, the power 
spectra for slow modes do not show the Kolmogorov slope. The slope is close to —2, 
which is suggestive of shock formation. At this moment, it is not clear whether or not 
the —2 slope is the true slope. In other words, the observed —2 slope might be due to 
the limited numerical resolution. Runs with higher numerical resolution should give the 
definite answer. 

Formation of shocks is expected as within slow modes for sufficiently high Mach num- 
bers the turbulence gets superAlfvenic. However, it is interesting that the turbulence 
anisotropy is definitely affected by magnetic field. 

Alfven modes follow the GS95 scaling as in the incompressible MHD turbulence. Alfven 
perturbations cascade to small scales over just one wave period, while the other non-linear 
interactions require more time. Therefore we expect that the non-linear interactions with 
other types of waves should affect Alfvenic cascade only marginally. Moreover, as the 
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Figure 2. M s ~ 10, Ma ~ 0.5, /3 ~ 0.01, and 512 3 grid points. Solid lines are for velocity 
spectra and dotted for magnetic field, (a) Spectra of Alfven modes follow a Kolmogorov-like 
power law. (b) Eddy shapes (contours of same second-order structure function, SF2) for velocity 
of Alfven modes shows anisotropy similar to the GS95. (rii oc r 2 / 3 or ku oc k 2 / 3 ). The structure 
functions are measured in directions perpendicular or parallel to the local mean magnetic field 
in real space. We obtain real-space velocity and magnetic fields by inverse Fourier transform 
of the projected fields, (c) Spectra of slow modes are a bit steeper than the Kolmogorov-like 
power law. They closer to k~ 2 . (d) Slow mode velocity shows anisotropy similar to the GS95. 
We obtain contours of equal SF2 directly in real space without going through the projection 
method, assuming slow mode velocity is nearly parallel to local mean magnetic field in low j3 
plasmas, (e) Spectra of fast modes show a possible departure from the Kolmogorov spectrum or 
the IK spectrum, (f) The magnetic SF2 of fast modes shows isotropy. 



Alfven waves are incompressible, the properties of the corresponding cascade do not 
depend on the sonic Mach number. 



4. Density statistics 

Density at low Mach numbers follow the GS95 scaling when the driving is incom- 
pressible (CL03). However, CL03 showed that this scaling substantially changes for high 
Mach numbers. Our high resolution results in Fig. 3 confirm the CL03 finding that at 
high Mach numbers density fluctuations get isotropic. Moreover, our present studies show 
that the spectrum of density gets substantially flatter than the GS95 one (see Fig. 4). 
Note, that a model of random shocks would produce a spectrum steeper than the GS95 
one. A possible origin of the flat spectrum is the superAlfvenic perturbations created 
by fast modes within density perturbations originated from slow modes. This particular 
regime is clearly identified in a review by CLV03 (see Fig. 9 therein). It may also be re- 
lated to the regime of superAlfvenic turbulence which is discussed in Norlund & Podoan 
(2003). However, alternative explanations of the shallow density fluctuations exist and 
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Figure 3. (a) Density structures for M s ~ 3. (b) Density structures for M. 

gets isotropic as the Mach number, M s , increases. 
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Figure 4. (a) Spectra for M s ~ 3 at t=3.7. (b) Spectra for M s ~ 10 at t=6.0. Dashed 
lines=density; solid=velocity; dotted=magnetic field. Density spectrum gets flatter as the Mach 
number, M s , increases. 



our ongoing work should clarify which process is actually responsible for the unusual 
density scaling that we observe. 

The flat spectrum of density perturbations in high Mach number turbulence is defi- 
nitely a very interesting phenomenon. It does affect the properties of interstellar medium, 
including the origin of the small scale clumps that are ubiquitous in the diffuse cold gas. 
This fiat density spectrum should affect the star formation processes. A study in Boldyrev 
et al. (2002) indeed testifies that the density spectrum is flat in molecular clouds. How- 
ever, a more systematic study is necessary. 



5. Discussion 

In the paper above we have studied scaling relations for Alfvcn, fast and slow modes. 
First question that comes to mind is how good is to speak about different modes in 
strong MHD turbulence. It is well known, for instance, that for strong perturbations 
various linear modes are coupled. Well, our work shows that the coupling is indeed a 
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very important phenomenon at the large scales where the total Mach number M tot , which 
normalizes the velocity perturbations by the sum of Alfvenic and sound velocities, is of the 
order of unity. At smaller scales when M to t becomes less than unity the perturbations that 
we identify through our separation of mode procedure develop without much interaction 
between fast and Alfven modes. This can be understood if we take into account that 
Alfvcn modes cascade just over one eddy turnover time. Therefore there is not enough 
time for the non-linear interaction between fast and Alfvcn modes to take place. At the 
same time, the shearing of slow modes by Alfvenic turbulence does take place over a 
turnover time. Therefore the slow modes mimic Alfven mode anisotropy. A new effect 
that we start observing at high resolution is steepening of the slow mode spectrum. This 
is indicative that while shearing the Alfven modes create shocks within the gas associated 
with slow modes. This is not so surprising as the Alfven velocity is smaller within the 
slow modes in low beta medium. 

We have not studied here a case of turbulence that is globally superAlfvenic. This is 
motivated by our belief that unless magnetic field is dynamically important the turbu- 
lence will behave like ordinary hydro turbulence. At the scale when the Alfven velocity 
gets of the order of the perturbation velocity, we should encounter turbulence that is 
studied here. 

Why would we care about those scalings? How wrong is it to use Kolmogorov scal- 
ings instead? Dynamics, chemistry and physics of molecular clouds (see Falgarone 1999) 
presents a complex of problems for which the exact scalings may be required to a differ- 
ent degree. If we talk about dynamics of interstellar dust or propagation of cosmic rays, 
one must account for the actual scalings and couplings of different modes (see reviews by 
Lazarian & Yan 2003, Lazarian et al. 2003). There are other problems, e.g. turbulent heat 
transport where the exact scaling of modes seems to be less important (Cho et al. 2003) . 

6. Conclusion 

In the brief, we have studied statistics of compressible MHD turbulence with high 
resolution numerical simulation for highly (M s ~ 10) supersonic low-/? case. We provided 
the decomposition of turbulence into Alfven, slow and fast modes. We have found that 
GS95 scaling is valid for Alfven modes: 

Alfvcn: E A (k) oc fc" 5/3 , k\\ oc k 2 / 3 . 

Slow modes also follows the GS95 anisotropy for highly supersonic low-/3 case: 

Slow: fcii oc k 2 / 3 . 

The velocity spectrum of slow modes tends to be steeper, which may be related to the 
formation of shocks. 

Density shows anisotropy for the mildly supersonic case. But, for the highly supersonic 
case, density shows more or less isotropic structures. We suspect that shock formation is 
responsible for the isotropization of density. Our present studies show that the spectrum 
of density gets substantially flatter than the GS95 one (see Fig. 4) . 
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